Mon. Not. R. Astron. Soc. , 000-000 (1996) 



A catalogue of galaxy cluster models 



Eelco van Kampen ' and Peter Katgert 

1 Sterrewacht Leiden, Postbus 9513, 2300 RA Leiden, The Netherlands 

2 Royal Observatory Edinburgh, Blackford Hill, Edinburgh EH9 3H.J 



Accepted ... Received ...; in original form ... 



ABSTRACT 

We present a technique to construct a fair sample of simulated galaxy clusters, and 
build such a sample for a specific cosmological structure formation scenario. Conven- 
tionally one extracts such a sample from a single low-resolution large-scale simulation. 
Here we simulate the clusters individually at high resolution. This is made possible 
by the method of constrained random fields, in which one can put linear constraints 
on peaks in the initial smoothed density field. 

We assume that these peaks are the progenitors of present-day rich clusters, and 
select clusters for a catalogue by selecting their initial peak parameters. We find 
that the final cluster mass can be well approximated by a linear function of both 
the amplitude and the curvature of the initial density peak. Because the probability 
distributions of these peak parameters are known, we can construct a model catalogue 
selected on expected final cluster mass. Such a catalogue will not have a well-defined 
richness limit, because the relation between richness and mass is fairly broad. However, 
by applying the appropriate completeness corrections, the results for the mass-selected 
catalogue can be compared with observations for richness-selected cluster catalogues. 

Each cluster model is evolved from its constrained initial conditions by means of 
an N-body integrator. This includes an algorithm for galaxy formation, so we produce 
two-component models consisting of dark matter background particles and galaxies. 
The latter allow us to directly obtain the observable properties of the cluster models, 
and match these to the observed properties to define the present time in the models, 
and thus derive the amplitude of the initial density fluctuation spectrum, erg. 

We build a model cluster catalogue for the Oo = 1 CDM scenario that is designed 
to mimic the EN A GS sample of rich Abell clusters. We use the distribution of richness, 
corrected for incompleteness, to fix the present epoch. We find <t 8 = 0.4 — 0.5, which 
is consistent with other determinations. The catalogue is 70 per cent complete for 
a richness larger than 50, but we do have a complete subsample for richnesses larger 
than 75. 

As a first test we compare the cumulative distribution of line-of-sight velocity 
dispersions to those found for several observational samples, and find that they match 
best for eg ~ 0.4. This means that we find consistent values for erg for the CDM 
flo = 1 scenario on cluster scales. 

Key words: galaxies: clustering - galaxies: formation - galaxies: evolution - cosmol- 
ogy: theory - dark matter - large-scale structure of Universe 



1 INTRODUCTION 

Cluster of galaxies are the largest bound structures in the 
Universe, but they are dynamically young and still contain 
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signatures of their initital conditions. Statistical proper- 
ties of samples of galaxy clusters can therefore provide con- 
straints for models of large-scale structure formation in the 
universe. Large catalogues of observed clusters of galaxies 
have been compiled, of which the Abell, Corwin & Olowin 
(1990, ACO from here on) catalogue is the best known, but 
complete X-ray selected cluster samples have recently be- 
come available (eg. Ebeling et al. 1996). The ACO clusters 
have been most extensively observed in the optical, as that 
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is where they have been defined. A renewed interest in their 
observational properties has recently led to a flood of new 
data, and to a better understanding of their properties. It 
is therefore desirable to construct a model catalogue that 
mimics observed ones, thus allowing a meaningful statisti- 
cal comparison between observations and predictions. In 
this paper we discuss a technique to construct such a model 
catalogue, and we build a catalogue for one particular cos- 
mological scenario, viz. f2o = 1 standard CDM. 

One approach for obtaining a sample of model clusters 
is to extract clusters from large-scale N-body simulations. 
This has been done by Frenk et al. (1990), and more re- 
cently by van Haarlem, Frenk & White (1996), and Eke, 
Cole & Frenk (1996). The advantage of this approach is 
that completeness of the resulting catalogue is automatically 
guaranteed. However, a disadvantage is that the resolution 
of the extracted cluster models is relatively poor since voids 
and filaments, which are simulated as well, take up most 
of the volume. This means that the fraction of particles in 
clusters is only about 10 per cent. We avoid this problem 
by running an ensemble of individual, high-resolution clus- 
ter simulations. This technique has been used by various 
authors in order to study the influence of the choice of ini- 
tial conditions on the final cluster properties (Evrard 1989; 
Evrard, Mohr, Fabricant & Geller 1993; van Haarlem & van 
de Weygaert 1993). However, because these authors did not 
aim to construct a fair sample, a statistical comparison with 
observed catalogues as performed by Frenk et al. (1990) and 
others was not feasible. This is exactly the goal we aim for 
here: to construct a statistically fair catalogue of individual 
cluster models, each of which is simulated at relatively high 
resolution. In this context the term 'fair sample' is meant to 
indicate a sample which is representative of a volume- limited 
sample. 

We add to this a method which includes the formation 
and merging of galaxies in dissipationless N-body simula- 
tions (van Kampen 1996). The basic idea is that each group 
of particles that is roughly in virial equilibrium is replaced 
by a single particle with mass, position, velocity and soften- 
ing corresponding to that group. This is done several times 
during the evolution. Since groups that form later in the 
evolution can contain one or more earlier formed galaxies, 
both growing and merging of galaxies are included in this 
simulation method, although somewhat crudely (improve- 
ment on this scheme is in progress). One reason for using 
this method is that groups of particles that can be identified 
as galaxies need to be protected from destruction within the 
main cluster by numerical two-body disruption (Carlberg 
1994, van Kampen 1995). But it also is highly desirable 
to be able to directly compare the properties of the mod- 
elled galaxy distribution to observations, instead of making 
assumptions on the relation between galaxies and dark mat- 
ter. 

Both the fluctuations in the galaxy distribution and in 
the underlying dark matter distribution grow with time. 
The relation between the two (usually described by some 
bias parameter), may change with time. In the simulations, 
the present epoch is identified as the time when the proper- 
ties of the simulated galaxy distribution match the observed 
ones. We use the normalisation obtained by van Kampen 
(1996) from a set of field models run for the same cosmology 
as adopted here. 



In order to construct a sample of individually simu- 
lated cluster models, we need to make a few assumptions. 
We assume that clusters of galaxies, as we observe them at 
the present epoch, formed from peaks in the initial density 
field smoothed at the length-scale appropriate for clusters. 
These initial peaks have several characteristics which can 
be used to predict whether they will evolve into rich Abell 
clusters. We also assume that the initial density fluctuations 
are Gaussian distributed. This allows the use of the theory 
of Gaussian random fields (Rice 1954; Doroshkevich 1970; 
Adler 1981; Peacock & Heavens 1985; Bardeen et al. 1986, 
BBKS from here on) which provides the probability distri- 
butions for the peak parameters in the early (linear) stages 
of the evolution of the density field. 

The construction of a sample of cluster models involves 
the generation of a set of initial conditions that produce a 
set of simulated clusters with the same statistical properties 
as for an observed sample. Such realisations of initial matter 
distributions that produce a cluster with specific properties 
can be obtained by means of the Hoffman-Ribak method of 
constructing constrained random fields (Hoffman & Ribak 
1991). We used the implementation of this method by van de 
Weygaert & Bertschinger (1996). An important limitation 
of this method is that it can only constrain linear function- 
al of a field. This means that the defining quantity of the 
catalogue has to be a linear functional as well. Given this 
limitation on the constrained random field method of gen- 
crating initial conditions, we argue (in Section 3.2) that the 
final cluster mass is the best choice for the defining quantity. 

In order to construct a catalogue, we assume that A/" rich 
Abell clusters within a given volume originate from those M 
initial density peaks that have the largest predicted final 
mass. We use the amplitude and the curvature of the ini- 
tial peak to predict the final cluster mass (Section 3.3 and 
Appendix B). This combination appears to give a better pre- 
diction of the actual final cluster masses than that based on 
the initial amplitude alone, which would be sufficient if the 
evolution would remain linear. Using the probability distri- 
butions from Gaussian random field theory we then sample 
the distribution of expected total masses to construct a fair 
statistical sample. 

In order to build a model catalogue, we need to choose 
a cosmological scenario, and an observational catalogue to 
compare to. We adopt the standard CDM 0,^ = 1 scenario, 
because it allows us to compare with other studies (Frenk 
et al. 1990; White, Efstathiou & Frenk 1993; Eke, Cole & 
Frenk 1996). It also has the property that the amplitude 
of the dark matter fluctuations need not be fixed at the 
outset. This is important, since we do not know a priori 
how the galaxy distribution will evolve with respect to the 
dark matter. 

The observational catalogue we try to mimic is a com- 
plete subset of the Abell catalogue for which data was gath- 
ered in the context of an ESO Key-programme (Katgert et 
al. 1996; Mazure et al. 1996; den Hartog 1995). This sample, 
denoted as ENACS from here on, is the largest well-defined 
catalogue presently available with extensive information on 
the internal structure of clusters of galaxies. 

Because there is bound to be significant scatter in the 
relation between richness and the expected final cluster 
mass, our model catalogue, constructed to be complete in 
mass, will not be complete in richness. In Section 5 we dis- 
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cuss the cause of this incompleteness, and apply a statistical 
method to correct for it. This gives us complete distribution 
functions for cluster properties like richness and line-of-sight 
velocity dispersion. 

In Section 6 we present the model counterpart of the 
ESO-Kpgm catalogue for the Qq — 1 CDM scenario. We use 
the distribution of cluster richness in addition to the galaxy 
auto-correlation function (van Kampen 1996) to fix the nor- 
malisation (in terms of as)- We find that both measures 
provide a consistent timing for the ilo — 1 CDM scenario. 

We can use the catalogue to test the validity of the 
chosen scenario by means of other measures, like the distri- 
bution of cluster shapes or line-of-sight velocity dispersions. 
We describe several such tests in Section 7. A preliminary 
comparison to observations and a more detailed study of the 
intrinsic properties of the model clusters can be found in van 
Kampen (1994). 



2 RICH ABELL CLUSTERS AND PEAKS 
IN GAUSSIAN RANDOM FIELDS 

Our aim is to build a sample of model clusters that can 
be compared to an observed catalogue. We therefore need 
to find a relation between the defining criteria for observed 
and model clusters. The former are found within the dis- 
crete distribution of galaxies as observed in projection on 
the sky, while the latter are defined in the continuous den- 
sity field which we sample with the particle distribution in a 
numerical simulation. In this section we discuss the obser- 
vational and theoretical ways to define clusters of galaxies, 
and also the numerical methods to build cluster models. 



2.1 Abell's observational definition of a rich 
cluster 

We will compare our models to rich Abell (1958) clusters, 
so we first list Abell's main criteria for the definition of an 
observed richness class > 1 cluster. They read: 

Richness criterion - 'A rich cluster must contain at least 
50 members that are not more than 2 magnitudes fainter 
than the third brightest member'. 

Compactness criterion - 'A rich cluster must be suffi- 
ciently compact that its 50 or more members are within a 
given (projected) radial distance R of its centre'. 

Abell's choice R = l.ShT 1 Mpc is usually referred to 
as the 'Abell radius'. The richness and compactness criteria 
jointly set the length- and mass-scale of Abell clusters. We 
now know that the Abell radius is somewhat small compared 
to more physical scale-lengths like the turnaround radius, so 
the Abell radius defines just the central region of a cluster. 
Most observations are restricted to an area within this ra- 
dius, or to an even smaller region (the cluster catalogue ex- 
tracted from e.g. the APM survey uses half the Abell radius, 
Dalton et al. 1992). 

Instead of a richness class we will use the richness mea- 
sure CacO) which actually gives the number of galaxies sat- 
isfying both of Abell's criteria. 



2.2 Clusters as peaks in the smoothed density 
field 

Theoretically, clusters are identified in the density contrast 
field 

P« - (P(r)> 



S(r) 



<p(r)> 



(1) 



as fluctuations at a certain mass scale (e.g. Kolb & Turner 
1990). This requires smoothing, or 'windowing', of the field 
at a corresponding length-scale. The smoothed density con- 
trast field is 



8 w (r,R w )=4n / «5(|r - r'|) W{r' , R w )rAr , 



(2) 



where W(r, Rw) is a normalised spherically symmetric win- 
dow function. Clusters are identified as the peaks in this 
smoothed density field, and the properties of these peaks 
characterise and define them. However, galaxy clusters 
exhibit highly non-linear evolution. This means that the 
present-day properties can not be linked directly to the ini- 
tial conditions from which they formed: we cannot simulate 
backwards in time. We are forced to make an educated guess 
for the initial conditions, simulate the formation of clusters 
starting from these conditions, and compare the resulting 
models to observations. This means that model clusters are 
necessarily defined in the initial density field, preferably at 
the epoch before clusters 'break away' or 'turn around' from 
the general expansion of the universe. 

To generate a catalogue of cluster models, we must de- 
fine models by peak parameters in the initial density field 
that can be sampled from a known probability distribution. 
If we choose peak parameters that are linear functionals of 
Gaussian random fields filtered with a Gaussian filter, we 
can use the results of BBKS. We therefore choose the Gaus- 
sian filter 



W G (r,R G ) = (2nR 2 G )-h- r2 ^ , 



(3) 



and make the implicit assumption that structure forms from 
small Gaussian density fluctuations in the early universe. 
Conveniently, integrals with a Gaussian kernel behave well 
and can often be found analytically. 

Note that the probability distributions derived by 
BBKS are only valid for the initial, linear regime (i.e. 
8\y <C 1). Non-linear evolution changes these distributions 
in an unknown manner, which makes the relation between 
the initial peak parameters and the final cluster model rather 
uncertain. Fortunately, the smoothing scale needed for clus- 
ter sized peaks is large enough that the evolution of clusters 
in the smoothed density field is not too non-linear, so that 
the BBKS distributions can be used. 

2.3 Initial peak parameters 

2. 3. 1 Definitions 

It is convenient to express peak parameters in terms of global 
parameters for the cosmological scenario in which a clus- 
ter is to be modelled. The general density fluctuation field 
smoothed with a Guassian with length-scale Rq can be char- 
acterised by the spectral moments 



2.i 



dk 
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where is the Fourier transform of the density field. Two 
special spectral parameters are usually defined: 



o-o(-Rg)o-2(-Rg) 



and R*(R G ) = V3- 



(5) 



They relate to the width of the spectrum and its correlation 
length respectively. For this paper we adopt the standard 
adiabatic CDM scenario for Slo = 1 with Rq = 4/i _1 Mpc, 
which has 7(4/i~ 1 ) = 0.729 and R*(4h~ 1 ) = 5.16/i _1 Mpc. 

We define the origin of the coordinate system at the 
centre of the peak, and define all peak parameters at this 
position. The primary parameter of a peak is its amplitude 
Sq(0, R G ), expressed in units of the r.m.s. fluctuation at the 
same scale: 



6 G (0,R G ) = ua (R G ) 



(6) 



We can further specify the peak using the first and second 
spatial derivatives of 5 G (r, Rq)- The gradient VS G (0,R G ) 
is zero for an extremum. This accounts for three parameters. 
The second derivative tensor ViVj8 G has six independent 
elements, resulting in six additional parameters. They char- 
acterise the shape (two parameters), the orientation (three 
parameters) and the size (one parameter) of the peak. The 
latter is the trace of ViVjS G , and is usually defined in a 
dimensionless form as the curvature: 



x(r,R G ) = - 



X7 2 8 G (r,R G ) 
°2(R G ) 



(7) 



This curvature (a measure for the width of the peak) is 
positive for a peak, and negative for a dip. The shape of 
the peak can be characterised by the two axial ratios of 
the isodensity contours around the peak, as found from the 
small r approximation of the smoothed density field around 
the peak using a second order Taylor expansion in spherical 
coordinates (r, 0, 4>) : 



r 2 2 

5 G (r) = S G (0) + yV 2 5 G (O)^ 5 (ai2,ai3,0,< 



(8) 



(BBKS). Here ^5(012, ^13, 0, <fr), given by equation (A7) in 
Appendix A, describes the asphericity, and 

- a - a in\ 

«12 = 7 , a i3 = - (9) 
c 

are the axial ratios, where a is the major axis, b the inter- 
mediate axis, and c the minor axis. 

2.3.2 Probability distributions 

BBKS give the probability distributions for the peak pa- 
rameters (the amplitude v, the curvature x, and the shape 
parameters a\i and 013) in the linear field smoothed with a 
Gaussian window of size R G . 

In order to make dependencies transparent, we use func- 
tions T (adapted from BBKS) which are listed for reference 
in Appendix A. The differential peak amplitude probability 
density is given in the form of the expected number of peaks 
within a volume V with an amplitude in the range [v, v + dv]: 



(10) 



The other probabilities are conditional; given a peak with 
height v, the probability distribution for its curvature to be 
equal to x is 



P{x\v,R G ) 



e T 2 (x) 
^2tt(1- 72)^1(7,^) 



(11) 



For a peak with curvature x, the axial ratios a\2 and a\z 
have a joint probability distribution 

P(a 12 ,a 13 \x) = gg£ ^3(ai2,ai3)x 8 e - §3; V^ 4 ( ai2 , ai 3) {u) 
0r r2\x) 

This distribution does not depend on v and R G directly, but 
only through x. All probabilities do depend on the spectral 
parameters j(R G ) and R*(R G ). 

2.4 Building individual cluster models 

2-4-1 Constrained initial conditions 

It is straightforward to generate initial configurations for 
an average patch of universe in a given cosmological sce- 
nario (e.g. Efstathiou et al. 1985). However, for our purpose 
we need initial conditions that will produce a cluster de- 
fined by certain peak parameters. This is made possible 
by the method of constrained random fields, pioneered by 
Bertschinger (1987). We use a code written by van de Wey- 
gaert & Bertschinger (1996), based on the Hoffman & Ribak 
(1991) method, which is very well suited for our purpose. It 
can sample an unsmoothed density field constrained to form 
a peak characterised by linear functionals of the Gaussian 
smoothed density field. Therefore the only restriction of this 
scheme is that we must use linear peak parameters to de- 
fine clusters in the initial conditions. This is an important 
limitation on the choice of a quantity that defines a model 
cluster catalogue (Section 3.2). 

2.4-2 Non-linear evolution and galaxy formation 

The highly non-linear evolution of clusters of galaxies re- 
quires the use of N-body methods. However, N-body sim- 
ulations suffer from a problem that is particularly severe 
for galaxy clusters: small groups of particles that represent 
galaxies in such simulations get disrupted by numerical and 
tidal effects (Carlberg 1994; van Kampen 1995). This prob- 
lem can be solved by replacing those groups of particles by 
a single 'galaxy particle' just after they have formed into a 
virialised system that resembles a galactic halo. This new 
galaxy particle is softened according to the size of the group 
it replaces, and gets the position and momentum of the 
centre-of-mass of the original group. In van Kampen (1995) 
such galaxy particles were only formed at a single epoch. 
Here an extension of this scheme to 'continuous' galaxy for- 
mation, which consists simply of applying the algorithm sev- 
eral times during the evolution, is used (van Kampen 1996). 
Since groups of particles can contain already formed galax- 
ies, merging is included as well, although somewhat crudely. 

The actual N-body code we use is the Barnes & Hut 
(1986) treecode, slightly adapted for the modelling of clus- 
ters (see van Kampen 1995 for details) and supplemented 
with the galaxy formation algorithm. For all simulations we 
generated initial conditions in a 64 3 cube with a size of 32 
/i _1 Mpc, using an implementation of the Hoffman & Ribak 
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(1991) method by van de Weygaert & Bertschinger (1996). 
Because we use non-periodic boundary conditions (the uni- 
verse outside the simulation volume is homogeneous), we cut 
a sphere out of the volume containing 64 3 particles. The ra- 
dius of this sphere is chosen such that it is sufficiently large 
to easily include the final turnaround radius, thus contain- 
ing the full cluster and its region of influence (see Appendix 
C for a discussion on the reliability of the simulations) . 

2.4-3 Model cluster richness 

Abell's richness measure Caco is defined in a magnitude 
interval. This means that we should correctly model the 
shape and the amplitude of the luminosity function in order 
to obtain a sensible richness measure from our numerical 
models. Frenk et al. (1988) found for CDM models that 
the shape of their mass function differs from that of the ob- 
served luminosity function. However, such results strongly 
depend on how one identifies galactic haloes at the present 
epoch (Suginohara & Suto 1992). The influence of the de- 
tails of the modelling of galaxy formation is also present in 
our simulations. We expect to get a more realistic lumi- 
nosity function since we identify and form galaxies during 
the simulation, rather than only at the end. If we assume 
a constant mass-to-light ratio, which should be reasonable 
for the factor of 30 in mass range, we find that the resulting 
luminosity function fits the observed one quite well. This 
means that we can obtain a realistic richness measure for 
our cluster models. 

The assumed constant M/L implies that the mass inter- 
val corresponding to Abell's richness criterion ranges from 
the mass of the third most massive galaxy to a mass 6.31 
times smaller. There is one property of our galaxies that 
plays a role here: the galaxy particle masses are multiples 
of the N-body particle mass, which is 7 x 10 10 M for our 
simulations. Before we determine the richness, we add or 
subtract a random fraction of half the N-body particle mass. 
This removes the artifical discreteness sufficiently well. 

Note that we do not have sufficient foreground and 
background galaxies in our simulations to mimic the ACO 
richness measure. For this reason we will use a different 
richness measure, which is statistically equivalent to Abell's 
(see Section 5.2). 

3 DEFINITION OF A MODEL 
CATALOGUE 

3.1 From initial density peaks to rich Abell 
clusters 

In order to build a catalogue of individual cluster models, we 
need to assume an one-to-one mapping of the peak parame- 
ters v and x in the initial smoothed density to properties of 
clusters in the final unsmoothed density field. However, this 
mapping might be complicated by initial peaks merging into 
a single final cluster, a worry that comes in two flavours. 

3.1.1 Mergers of rich cluster peaks 

In a hierarchical structure formation scenario, like CDM, 
merging of clumps into larger clumps is an ongoing process. 
Because we assume a one-to-one mapping of peaks in the 



initial density field to overdensities in the final density dis- 
tribution, we need to check if merging is still important on 
cluster scales. If evolution would remain linear, peaks that 
were separate in the initial smoothed density field would re- 
main separate for the same comoving smoothing scale. They 
can only merge due to non-linear evolution of the density 
field. We estimate the fraction of clusters that can merge 
from the statistics of the smoothed density field, i.e. the cor- 
relation function of Abell cluster peaks, and their peculiar 
velocities. The average peculiar velocity for a cluster-sized 
peak is about 300 km s _1 , with a maximum of about 600 
km s _1 , for a D.q = 1 CDM cosmology with o% = 0.6 (de 
Theije, van Kampen and Slijkhuis, 1996). In this scenario, 
a typical cluster is therefore expected to travel about 2h~ 1 
Mpc up to the present epoch. 

The average nearest-neighbour distance r nn is given by 

o /Tnn 

rnn = -r < nc > _1 -3 / £cc(r)r 2 dr (13) 

47r Jo 
where n c is the cluster number density, and £ C c(r) the 
cluster-cluster correlation function, which is observed to be 
a power-law, 

£cc = (r/rcc) -1 - 8 , (14) 
where the cluster-cluster correlation length r cc is around 
I8/1 -1 Mpc (eg. Peebles 1993). This means that for a mean 
cluster density of 8.6 x 10~ 6 /i 3 Mpc~ 3 (Mazure et al. 1996) 
the average nearest-neighbour distance is 20h~ 1 . Therefore 
one does not expect two Abell sized peaks to merge within 
a Hubble time. This is confirmed from an estimate of the 
probability of finding, at the present epoch, a cluster-cluster 
separation of less than 2fe _1 Mpc (the typical distance trav- 
elled up to the present epoch), using the observed £ cc . We 
derive a probability of 2 per cent, i.e. about 2 clusters in the 
catalogue. This means that we expect one close pair that 
could merge if moving towards each other. Note that we 
derived this probability for the present epoch, and that it 
must have been smaller at earlier times. 

3.1.2 Late cluster formation from subcluster mergers 

A cluster can also form from two or more subclusters, which 
are more abundant than rich clusters themselves, through 
a merger at late times. One needs somewhat special initial 
conditions for clusters to form this way, because the separa- 
tion of such subclusters in the initial conditions cannot be 
too large or they would remain distinct, despite non- linear 
clustering. It is more likely that late merging involves un- 
equal mass clumps, where the main clump would not have 
been massive enough to become a rich cluster were it not 
for the 'secondary infall' of nearby subclumps. 

Such cases of cluster formation are unavoidable in hier- 
archical formation scenarios. The direct way to access their 
significance is to look at a large-scale simulation of volume 
equal to that of the whole catalogue (see eg. van Haarlem, 
Frenk & White 1996). We will discuss this question in Sec- 
tion 5. 

3.2 The catalogue defining quantity 

3.2.1 Total cluster mass 

As a single cluster model is defined by its local initial condi- 
tions, a model catalogue of clusters is defined by the statis- 
tics of peaks in the initial cosmological density field. We 
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will construct a model catalogue under the assumption that 
a certain subset of peaks in the initial smoothed density field 
will each produce a cluster at the present epoch. In order to 
build a catalogue that mimics a sample of rich Abell clusters, 
we need to establish a relation between peak characteristics 
in the initial density field and the criteria used by Abell at 
the present epoch. Unfortunately, his criteria are not eas- 
ily translated into the formalism of Gaussian random fields. 
Therefore we are left with a choice for a sensible physical 
or observable quantity that does provide a relation between 
initial and final conditions. 

Physical parameters are very easily obtained for model 
clusters since the data are clean and unprojected. Far less 
choice is available in the observations, which provide pro- 
jected, biased and generally less clean data. We will fo- 
cus mainly on galaxies, not on dark matter, although it is 
probably the dominant mass component. The dark matter 
distribution can only be obtained through indirect measure- 
ments like the X-ray brightness distribution (e.g. Sarazin 
1986; Forman & Jones 1990) or the effects of gravitational 
lensing (e.g. Fort & Mellier 1994). 

Which physically significant quantities can we derive 
from the galaxy distribution? The quantities that are best 
observable are the two spatial coordinates on the sky, which 
can be obtained very accurately. Nevertheless, the position 
of the centre of the projected galaxy distribution, required 
for most of the derived properties, is relatively difficult to 
obtain. Global properties derived from the projected dis- 
tribution are the total number of galaxies within a certain 
radius, which was used by Abell for his richness criterion, 
the shape (ellipticity) and orientation of the distribution 
(e.g. Binggeli 1982; Rhee, van Haarlem & Katgert 1989; 
de Theije, Katgert & van Kampen 1995), measures of sub- 
structure (reviewed by Beers 1992), the azimuthally aver- 
aged surface density profile (e.g. Beers & Tonry 1986; Rhee 
et al., ibid), and the projected mean harmonic radius, which 
is often used for mass estimates. 

With regard to velocity information, only line-of-sight 
components are readily observable. Apart from the distance 
to the cluster, the most useful global quantity derived from 
these is the central line-of-sight velocity dispersion. This 
would be a good candidate for the catalogue defining quan- 
tity, except that it is a non-linear quantity, and can therefore 
not be used in the the Hoffman- Ribak method of constrained 
random fields. The total velocity dispersion is correlated 
with the total mass of the cluster, which is also one of its 
most basic physical properties. We demonstrate below that 
the final cluster mass can be estimated fairly reliably from 
the initial peak parameters, and we therefore choose the ex- 
pected final mass to be the defining quantity of a catalogue 
of cluster models. 

Although in principle mass estimates can be obtained 
from the combination of the projected galaxy distribution 
and the central line-of-sight velocity dispersion, reliable total 
masses are not known for most Abell clusters. All we can 
do is to hope that richness correlates sufficiently well with 
mass, and that those peaks in the Gaussian smoothed initial 
density field that have the largest expected final mass are 
the progenitors of the rich Abell clusters. 



3.2.2 Smoothing scales 

Before we discuss the relation between the final mass of a 
cluster and the parameters of the initial density peak, wc 
need to choose the smoothing scale that is appropriate for 
clusters. We base our choice on an argument involving the 
cluster mass. In general, the total mass for a given spherical 
window function W(r, Rw) i s given by 

M w = V w p[l + S w (0,R w )], (15) 

where the volume Vyy(Rw) of the window function is 

V W (Rw) =4tt J W(r,R w )r 2 dr . (16) 

The (observational) compactness criterion implies that 
Abell clusters are defined using a Top-Hat window 

with i?TH = l-5h~ Mpc, the Abell radius. The mass en- 
closed by a sphere this size is about 4 x 10 h~ Mq for the 
background universe (for Qo — 1). A 'typical' cluster has a 
mass of about 3 x 10 14 n ft _1 Mo (e.g. Peebles 1993), so the 
cluster mass originates from a volume 70 times larger. This 
corresponds to a Top-Hat filter of about 6ft _1 Mpc in the 
initial density field. The corresponding Gaussian smoothing 
radius Rq is found by setting the enclosed volumes equal 
for both windows (BBKS). This gives Rq w 4/i~ 1 Mpc, as 
the ratio between the two radii is (4/187T) 1 / 6 w 2/3. 

3.3 Cluster mass threshold 

In this section we estimate the total mass of a cluster given 
the properties of the peak in the initial smoothed density 
field it originates from. 

3.3.1 Peak amplitude threshold: 'linear' mass 

In the linear approximation the final total mass is propor- 
tional to the peak amplitude in the initial smoothed density 
field. We define model clusters using the Gaussian window 
function with radius Rq = 4/i _1 Mpc. If we take the same 
value for the present density field as for the initial field, and 
assume the field to evolve linearly, the total mass is deter- 
mined exactly by the peak amplitude v alone and given by 

Mq{v) = (2nR 2 G f 2 p[l + vc-o] , (18) 

because the relation 8q — voq remains valid throughout the 
evolution. In that case the linear relation ao{t) = aao(to) 
describes the time dependence, where io is the present epoch 
and a is the cosmological expansion factor (with a (to) = 1), 
and v remains constant. 

We could then define a model catalogue by setting a 
threshold on v, which is that value of v for which iV pcak (> 
v ,Rq) = A/", where M is again the total number of clusters 
in the catalogue with volume V. The threshold f m ; n , which 
is formally an expectation value, should then correspond to 
the threshold in the defining quantity of the observational 
catalogue. In the case of rich Abell cluster catalogues this 
is the value of 50 for the Abell richness measure. 

We can calculate the expected maximum peak ampli- 
tude i/max in that volume from iV peak (> v,Rq) = 1. We 



A catalogue of cluster models 7 



10000.0 



1000.0 



o 
Q. 



100.0 



~ 10.0 



1.0 



0.1 





















CDM 








COM.!! 0.2 








k" 2 

CDM, h=l 




. \\ 





2 3 
v 



Figure 1. Expected range in v for 100 clusters in a volume of 
10 7 /i~ 3 Mpc 3 . All scenarios are for Qq = 1 and h = ^ unless 
stated otherwise. 



illustrate this in Fig. 1 for several cosmological parameters 
and fluctuation spectra, smoothed on the cluster scale of 
4ft~ 1 Mpc. The heavy horizontal line determines the thresh- 
old i/ min for the choices N = 100 and V = 10 7 /i~ 3 Mpc 3 , 
which are typcial values for a cluster catalogue. The thin 
line indicates the expected which is supposed to rep- 

resent the largest cluster in V . Of course, larger peak am- 
plitudes will occur in larger volumes. 

Note that it is not sufficient to just threshold the 'lin- 
ear' mass in the initial smoothed density field in order to 
guarantee the presence of a rich Abell cluster in the final 
density field. Because richness and mass are not likely to 
be perfectly correlated, there are peaks with a small 'linear' 
mass which produce rich clusters (and vice versa), making 
the catalogue selected on mass incomplete in richness. Since 
we cannot use the constrained random field method to con- 
strain Abell richness, we have to tighten the correlation of 
initial peak parameters and final richness as much as possi- 
ble. One way of doing that is to go from 'linear' mass, as 
given by the peak parameter v, to a 'non-linear' mass. 

Several non-linear effects come into play which deter- 
mine the final mass of a galaxy cluster. The collapse time 
of a cluster not only depends on the peak amplitude of the 
initial overdensity, but also on its curvature x (van Haarlem 
& van de Weygaert 1993), on tidal effects from the environ- 
ment, on the presence of substructure (Cavaliere et al. 1986) 
which can result in pre-virialisation (Peebles 1990) and thus 
a slower collapse, and on the amount of shear in the initial 
velocity field (Bertschinger & Jain 1994). If the non- linear 
evolution of these properties is significant, the peak in the 
final density field may not coincide with the peak in the 
initial density field, not even for the field smoothed on the 
scale of 4/i~ 1 Mpc. 

We show below that the effect of the peak parameter 
x, which is related to the slope of the initial density profile, 
can be taken into account, and can be used for a better 
approximation of the expected final cluster mass. 



3.3.2 Expected 'non-linear' cluster mass 

Cluster peaks in the density field smoothed at Ah~ l become 
non-linear. We find that Sq can grow up to about 10 for 
the structure formation scenario that we have adopted, viz. 
standard CDM for Q,q = 1. Therefore v becomes a function 
of time, and ao(t), the r.m.s. fluctuation on scale Rq, will 
become quasi-nonlinear. Note that u(t) is a locally evolv- 
ing quantity whereas ao(t) evolves globally, so the latter will 
grow less rapidly than the former which describes the evo- 
lution of an individual peak. Furthermore, the effect of the 
initial curvature of the peak, as set by x, will become im- 
portant. Van Haarlem & van de Weygaert (1993) found 
that the final cluster mass is larger for smaller x, i.e. more 
extended initial configurations. We can try to incorporate 
these effects by looking at the possible evolution of the ini- 
tial smoothed radial density constrast profile around a peak 
with given v and x. The initial density profile is given by 
BBKS: 

- a:/7 



S G (r) = 



v — 7a; 
«ro(l -7 2 ) 



&}(»■) + 



3a (l- 7 2 ) 



(19) 



where £c( r ) i s the smoothed density autocorrelation func- 
tion. Remember that all quantities depend on the smooth- 
ing radius Rq. In the linear regime Sq(0) = voq, and 
£g(0) — °o by definition, so we can derive that initially 



v'Cg(o) 



o 2 2 
-37 °~0 



(20) 



For small r both terms in eq. (19) initially contribute roughly 
equally to the density profile for most cosmological scenar- 
ios. For large r the V 2 £g( t term quickly vanishes and the 
profile is well approximated by just the first term. For small 
r the profile is well approximated by a second order Taylor 
expansion around the peak, as in equation (8). Here we will 
just consider the spherically averaged profile, which means 
that the contribution from J-5 cancels. 

We now substitute non-linear quantities in eq. (19) in 
order to get an evolution equation for the (smoothed) den- 
sity profile. If the evolution would remain linear, v, x and 
7 would remain constant, and o"o would evolve linearly with 
the cosmological expansion factor a. The autocorrelation 
function scales with <tq, and so do its derivatives because the 
slope of (, G (r,t) remains constant, so that the a;(i)-terms in 
eq. (19) cancel. How does non-linear evolution change this 
behaviour ? The r.m.s. fluctuation gq at 4/i _1 Mpc Gaus- 
sian smoothing is about unity at the present epoch for most 
cosmological scenarios, and its time evolution therefore be- 
comes non-linear. Yet, v and x become non-linear more 
quickly because they are defined for peaks. Now, if £ G (0, t) 
would grow at the same rate as V 2 £g(0, t), both a;(t)-terms 
in eq. (19) still cancel and we would only see non-linear be- 
haviour in <5g(0, t) due to i/(t)ao{t). 

However, non-linear evolution causes the slope of 
CgW to become steeper, which changes the evolution of 
V 2 £g(0, t) with respect to £g(0, t). So the second term of 
eq. (19) evolves somewhat faster than the first because of the 
increasing slope of £g(*)- This might seem unimportant, but 
the x-terms, which grow relatively fast, do not cancel any- 
more. This is demonstrated in more detail in Appendix B. 
If we use a first order approximation for v{t) and x(t), we 
find a relation of the form 

M G ( Vi ,x h t ) = (27ri? G ) 3/2 p[l + ( C o + Cl ^ + C2a;i )ao] ,(21) 
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where the cj are constants at epoch to, and v\ and x\ are 
the initial peak parameters. If the field and its derivatives 
would just evolve linearly, cq = 0, c\ = 1 and C2 = 0. Non- 
linear evolution results in a non- vanishing Co , where the sign 
depends on the exact balance of the constants in the linear 
approximations for u(t) and x(t), c\ > 1 and C2 < 0. Again, 
see Appendix B for more details. 

3.4 Sampling constraints above the expected 
final mass threshold 

We have defined a model catalogue of Af clusters as follows. 
We first draw a large number of (v, x) sets for a sufficiently 
small ^ m i n found from (10), and then apply the M(y,x) 
threshold so that N models remain. However, the values 
for the parameters Cj in (21) are unknown a priori because 
no model has actually been built yet. So one first needs 
to choose a few representative models from the N models 
selected for small deviations of the linear values cq = 0, c\ = 
1 and C2 = 0, for example Co = 0, c\ = 1.5 and C2 = —0.5 
(see Appendix B). After simulating these models, one can 
fit the Cj to the measured total mass in these simulations 
within a comoving sphere of radius 6/i _1 Mpc around the 
centre. This provides a new total mass threshold, which 
can then be used to adjust the selection of the M models 
from the large ensemble. If the initial sets of constraints 
are chosen sensibly, along the lines described in Appendix 
B, and not too close to the initial threshold, most of the 
trial models will remain in the catalogue. This iterative 
procedure is repeated until all J\f models are run and the Cj 
have converged to values that define the final catalogue. 

Each model is defined not only by v and x, but also 
by its shape parameters. These do not form part of the 
predicted mass relation that defines the catalogue, but arc 
still chosen so that we can also study the dependence of 
the final cluster models on initial shape. Given v and x, 
the parameters ayi and 013 are drawn simultaneously from 
the probability distribution in eq. (12). The orientation is 
drawn from an isotropic distribution, as the phases of the 
Gaussian density field are random. The actual orientation 
is important only to determine the projected shape of the 
object on the sky. 

Deviates are generated using the rejection method (e.g. 
Press et al. 1988). To draw v, a convenient comparison 
function is found by substituting v 3 for ^1(7,1/) in (10). 
The comparison function needed for drawing x is obtained 
by replacing ^{x) in (11) with x 3 (see also Appendix A). 
For the joint sampling of the shape parameters the compar- 
ison function is simply a constant that is larger than the 
maximum of (12). Deviates of a\i and 013 are (uniformly) 
generated between 1 (by definition) and some maximum for 
which the probability is really small already, for which we 
took 10. This may not be the most efficient algorithm, but 
it works sufficiently well for our purposes given that (12) is 
rather complicated. 

4 CONSTRUCTION OF A MODEL 
CATALOGUE 

4.1 Choice of the cosmological scenario 

We build a catalogue of model clusters within the standard 
CDM scenario, i.e. we adopt the fio = 1 adiabatic CDM ini- 
tial fluctuation spectrum (Davis et al. 1985). This scenario 
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Figure 2. Selected csomological density fluctuation spectra, 
plotted in the form of u p (k) = (47rfc 3 ) 1 / 2 |<5 fc |. 

has the advantage that the amplitude of the initial den- 
sity fluctuation spectrum, denoted by as, can be fixed after 
having run the models, by comparing the properties of the 
galaxy distribution to observational data. Fig. 2 shows the 
unbiased standard firj=l adiabatic CDM spectrum as given 
by Davis et al. (1985) in the form of a p (k) = (47rfc 3 ) 1/2 |<5 fe |, 
which is the r.m.s. contribution to the density fluctuation 
field. If a s is less than the observed value of unity, this de- 
creases the amplitude of the spectrum by a factor of a$ . for 
all wavenumbers. The vertical lines indicate the range of 
wavenumbers included in our simulations. The dashed line 
shows the Qq — 1 CDM scenario, which we adopted here, 
for g = 0.5, a value that should not be far off from the one 
that best fits the observational data. We also show a few 
other scenarios, with values for <jg such that the spectrum is 
not too different from the fio — 1 CDM spectrum within the 
range of wavenumbers present in the simulations. Of course 
they will be significantly different at larger and/or smaller 
wavenumbers: eg. the featureless power-law spectrum with 
index —2 (dashed-dotted line) has substantially more power 
on both larger and smaller scales. 

Another important reason for choosing the Q^t — 1 
CDM scenario is that it allows comparison to published re- 
sults, notably those of Frenk et al. (1990) and White et al. 
(1993). Furthermore, CDM does not have many more prob- 
lems than most of the competing scenarios (e.g. Davis et al. 
1993; Ostriker 1993; Bertschinger 1993), and explains many 
observations, although not in its original unbiased form. 

4.2 The target observational catalogue 

As we want our sample of model clusters to be fair, i.e. 
statistically complete for a specific volume, its observational 
counterpart also needs to be as complete as possible. A 
sample that was constructed with this requirement in mind 
is one for which data has been obtained in the context of an 
ESO Key-programme. This sample of rich Abell clusters, 
the ENACS sample for short, is based on a subsample of 
128 R > 1 clusters selected from both the Abell and ACO 
catalogues in a cone of about 2.6 steradian out to z = 0.1 
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Figure 3. Selection of models on expected total mass within 
4h~ ^Mpc at erg = 0.46 from 200 generated (y, x) sets. The thick 
solid line selects the 99 cluster models that have been simulated 
(filled symbols) and form the model catalogue. The dashed line 
denotes the expectation value of x for a given v (cq. AW). The 
thin solid lines show probability density contours of P(x\v) (see 
cq. (11) in Section 2.3, and Section 4.4.3). 



(Katgert et al. 1996). It occupies a single closed volume that 
is well away from the Zone of Avoidance. 

We choose a z < 0.08 subsample of the ENACS sample 
(which goes out to z = 0.10) as the observational catalogue 
to be mimicked, because it is the most complete subsam- 
ple selectable which is large enough for statistical purposes 
(Mazure et al. 1996). It contains 76 clusters in a volume 
of 9.2 x 10 6 /i _3 Mpc 3 . After corrections for completeness, 
Mazure et al. (ibid.) found that there should have been 
79 in that volume, implying a comoving number density of 
8.6 ± 0.6 x 10- 6 /i 3 Mp C - 3 for R > 1 Abell clusters. 




4.3 The Q = 1 CDM model catalogue 

We now define a catalogue of cluster models that mimics the 
ENACS sample. For the volume of the catalogue we take 
10 7 /i _3 Mpc 3 . This means that according to the number 
density found by Mazure et al. (1996) there should be 86 
rich Abell clusters in this volume, of which 83 define an 
ACO-like observational sample. To get close to the target 
of 83 cluster models, we built a simulation set that consists 
of 99 cluster models. This redundancy somewhat alleviates 
completeness problems associated with thresholding on the 
predicted final mass, which arises from scatter in the relation 
of richness versus final mass (Sections 5 and 6). 

We generated the simulation set for the CDM f2o = 1 
scenario using the technique described in Section 3.4. We 
used 200 sets of (v, x) from which the desired 99 will be se- 
lected. The 200 sets are shown in Fig. 3 in the form of a 
scatter plot of v versus x. The solid line depicts the fitted 
M(u, x) threshold found after the last iteration that sepa- 
rates the 200 sets in 99 candidate catalogue entries (filled 
circles), and 101 models that are below the expected final 
mass threshold (open circles). 



Figure 4. Fits for the expected masses Mq(u) (top panel) 
and Mq(v,x) (bottom panel) to the total mass found within 
6/i -1 Mpc in the simulations for erg = 0.46. 

The present epoch for the adopted scenario, as ob- 
tained from comparison to the galaxy autocorrelation func- 
tion, is most likely to be erg = 0.4 — 0.5 (Davis et al. 
1985, van Kampen 1994). Note that the COBE measure- 
ments require a rather different range of values, as is dis- 
cussed in Section 8. Our simulations have data saved for 
cr 8 = {0.25,0.31,0.36,0.41,0.46,0.51,0.56,0.61,...}, so we 
took <rg = 0.46 for the purpose of fitting the Cj in eq. (21). 
In Fig. 4 we show two fits: one for the linear mass estimate 
based on a linear function of v alone, and one for Mq(v,x), 
given by (21). The first clearly has a larger spread than 
Mq(u,x); the variance of the residuals of the fit, cr% s , is 
twice as large. We found cq — —4.3, c\ — 4.9 ± 0.3 and 
C2 = —1.3 ± 0.2, with no correlation between the residual 
scatter in the Mq (u, x) relation and the shape parameters 
a\2 and 013. We suspect that the remaining scatter is mostly 
due to the non-linear evolution and the difference between 
Gaussian and Top-Hat smoothing, as well as to differences 
in the initial velocity field, notably the shear (Bertschinger 
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& Jain 1994; van de Weygaert & Babul 1994). This latter 
effect will be investigated in another paper. 

4.4 Testing the construction technique 

With the aid of a conventional large-scale, low-resolution 
simulation we have tested the assumptions made in the 
method we use to construct a catalogue of galaxy cluster 
models. What we test here is whether the final galaxy clus- 
ters (i.e. the final peaks) indeed have a one-to-one correspon- 
dence to initial density peaks, and how well the expected 
mass estimator Mq performs. 

4-4-1 Description of the large-scale simulation 

We evolved a patch of universe in a cubic volume of 
(256/i _1 Mpc) 3 using a P M code (Bertschinger & Gelb 
1991), and saved particle positions and velocities at the same 
epochs as for the individual cluster models. A more detailed 
description of this simulation can be found in de Theije, van 
Kampen & Slijkhuis (1996). Given a cluster density for rich 
Abell clusters of 8.6 x lCT 6 /i 3 Mpc~ 3 , we expect 144 of such 
clusters in the chosen volume. We employed 128 3 particles, 
so that each particle has a mass of 4.4 x 10 12 M Q . This 
provides sufficient resolution to identify rich clusters and 
trace them back to the initial density field. Note that the 
resolution is significantly lower than for the individual sim- 
ulations, so that properties obtained for the clusters in the 
single large-scale simulation will be much noisier than for 
their individually simulated counterparts. 

After running the simulation, the 4/i -1 Mpc Gaussian 
smoothed density field and its second derivative where ob- 
tained on a 128 3 lattice for each output time. Peaks were 
identified on these lattices by searching for those cells which 
have a larger density than their 26 neighbouring cells. We 
then use second order polynomial fitting to find the off-grid 
positions of the peaks, and the actual peak parameters v 
and x at those positions. At the final epoch these peaks 
correspond to N-body particle groups which should be the 
final-epoch rich galaxy clusters. 

4-4-2 Tracing cluster peaks 

We first look at the link between initial and final density 
peaks by linking the former to the latter using our single 
large-scale simulation. We link the peaks at different epochs 
by finding the smallest spatial separation between the re- 
spective positions in consecutive time-slices. 

We found that none of the initial peaks have merged 
during the evolution up to erg — 0.6, well above the value of 
o g = 0.4 — 0.5 found previously (see Section 4.3). This justi- 
fies our use of initial peaks to construct a cluster catalogue, 
and confirms the estimate in Section 3.1.1. 

4-4-3 Distribution and evolution of v and x 

We examined the (y, x) distribution found in the single large- 
scale simulation, which is complete by definition, in order to 
test the generation of deviates for the peak parameters that 
define the catalogue of high-resolution models. We indeed 
find a similar distribution for (y, x) as for the models making 
up the model catalogue (see Fig. 3), so the sampling of the 
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Figure 5. (a) Relations between linear and final v (closed 
symbols) and x (open symbols) from the model catalogue, for 
erg = 0.46. (b) same as (a), but for those (x[,Xq) that have 
\v\ — 3 < 0.05. The dashed line is a fit to these points. Open tri- 
angles are for simulations by van Haarlem & van de Weygaert 
(1993), for erg = 1.0. A fit to these three points, which all 
have v\ = 3, is shown as a dotted line. Notice the evolution 
from solid line (erg = 0), via dashed line (erg = 0.46), to dotted 
line (eg = 1.0). (c) uz e \ and xz e \ versus the ratios fZel/^i an d 
x Zel/ x i)> with solid lines showing fits to these, which were used 
for Figs. 6(b) and 7. See text for more details. 
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Figure 6. (a) Evolution of the non-linear to linear ratio for the 
peak amplitude, u, and the peak curvature, x, for all catalogue 
models, (b) Same as (a), but for the peaks found in the single 
large-scale simulation. 

latter is representative (see also van Kampen 1994). 

For the catalogue of high-resolution models we show, 
in Fig. 5(a), for each simulation, the relation between the 
initial peak parameters, v\ and x\, and their final values vq 
and xrj a t the present epoch, to- Filled symbols are for v, 
whereas open ones are for x. The solid line indicates linear 
evolution, for which the peak parameters remain constant 
(i.e. vq — v\ and xq = x\). Whereas the evolution of v 
is mildly non- linear, the curvature x evolves well into the 
non-linear regime, but with a large scatter. 

In Fig. 5(b) we show the results of three simulations 
by van Haarlem & van de Weygaert (1993) (open triangles), 
who give initial-final relations for x for peaks with V[ = 3 
(note that they have their final epoch at as = 1, whereas we 
have as = 0.46). The dotted line is a fit to those points. We 
also include those models in Fig. 5(a) with \v\ — 3| < 0.05, 
as those are closest to the v\ = 3 peaks of van Haarlem & 
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Figure 7. Fit for the expected mass Mq of peaks in the single 
large-scale simulation to the total mass found within 6h~ 1 Mpc 
around these peaks. 

van de Weygaert (ibid.). A fit to their data is shown as a 
dashed line, and it can be seen that the slope of the initial- 
final relationship for x at a fixed v\ evolves from unity for 
the initial epoch (solid line), through zero near og w 0.4, to 
a steepening negative slope for larger ag. 

The evolution of the peak parameters of the high- 
resolution models might not be similar to that in the low- 
resolution simulation if the size of the high-resolution sim- 
ulation spheres was not chosen large enough, as that would 
artificially reduce large-scale tidal fields. We show the evo- 
lution with as of the ratio of final to initial values of v and 
x for each high-resolution cluster model in Fig. 6(a). 

A similar plot for the low-resolution simulation cannot 
be readily made, as we do not know the initial values of the 
peak parameters, but only v Zc \ and Xi e \, the values at the 
end of the Zeldovich' approximation which was used to cal- 
culate the linear evolution and the first part of the non-linear 
evolution. In Fig. 5(c) we show the relation between the 
value of the peak parameters at the start of the N-body in- 
tegration (vzd an d Kzeb an d the ratios ^Zel/^i and KZel/^i)- 
This figure shows that at the start of the N-body integra- 
tion the amplitude v of the peak is on, average, 1.1 times 
larger than the initial value, while the curvature is on aver- 
age 1.4 times larger than its initial counterpart. Using Fig. 
5(c), we have produced Fig. 6(b), which shows the evolution 
of the peak parameters for the low-resolution simulation. 
The initial scatter in Fig. 6(b) was put in by hand; it re- 
flects the scatter in the fits from Fig. 5(c) and was sampled 
from a Gaussian distribution. It appears that the individual 
simulations perform well; the peak parameters from the sin- 
gle large-scale simulation show a very similar evolution with 
time as those for the collection of high-resolution models, 
although they are somewhat noisier. 



4-4-4 Fairness of the expected mass relation 

Using the centre-of-mass of the actual N-body particle group 
near the density peaks in the low-resolution simulation, we 
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obtain the mass M s j m within a radius of 6.0 h~ Mpc (i.e. the 
Top-Hat mass-scale that was used to define the catalogue). 
We expect 144 rich Abell clusters in this volume (see above), 
so we select a threshold for Mq that produces that number. 

To test the Mq — M s ; m catalogue defining relation, we 
trace back clusters found at the present epoch to peaks in 
the smoothed initial density field. In Fig. 7 we plotted the 
final mass M s j m against the predicted final mass Mq that is 
calculated from the initial peak parameters according to eq. 
(21), where we again have used the fits of Fig. 5(b) to get 
from the 'Zeldovich' peak parameters to the initial ones. 

The relation in Fig. 7 compares quite well to the relation 
for the single, high-resolution, cluster models, as shown in 
Fig. 4(b). Although the constants cj are slightly different, 
the ratio C1/C2, which determines the slope of the final mass 
threshold, is almost the same (3.54 as compared to 3.76). 

5 COMPLETENESS 

5.1 Completeness of the ENACS sample 

The ENACS sample was constructed to be complete in rich- 
ness (Katgert et al. 1996) by selecting a subsample from the 
ACO catalogue within a particular volume (see Section 4.3). 
However, from comparison to the Edinburgh-Durham Clus- 
ter Catalogue (EDCC from here on), Mazure et al. (1996) 
estimated that the full R > 1 (Caco > 50) ACO catalogue 
is 94 per cent complete out to z — 0.08. We have to ac- 
cept this incompleteness, since the ACO catalogue is still 
the largest observational sample available. In the follow- 
ing we therefore assess completeness with respect to the full 
ACO cluster catalogue, keeping in mind that this catalogue 
is only 94 per cent complete. 

We see from Fig. 8 that the ENACS sample, as a sub- 
sample of the ACO catalogue, is 100 per cent complete out to 
z — 0.08. Plotted is the distribution of richness Caco °f au 
ACO clusters (thin solid line), ACO clusters out to z = 0.08 
(thick solid line) , and ENACS clusters out to the same red- 
shift (dotted line), where the latter is rescaled to the ACO 
volume. The distribution of the latter two is clearly statis- 
tically the same, indicating that the ENACS catalogue is a 
representative subsample of the ACO catalogue. An expo- 
nential function of the form n(R) = a e~ bR has been fitted 
to the richness distribution of the z < 0.08 ACO clusters for 
use in Section 6, and is shown in Fig. 8 as a dashed line. 

5.1.1 ACO and 3D richness 

Both Abell (1958) and ACO (1990) have corrected their rich- 
ness determinations for foreground and background galaxies, 
by subtracting a 'field correction' term from the raw number 
count. Since we model clusters individually, our simulation 
volume does not contain any obvious foreground or back- 
ground galaxies. Common observational practice is to use 
a 'gapper' criterion to remove such galaxies, i.e. eliminate 
all galaxies that are on the non-cluster side of a gap in the 
line-of-sight velocity distribution for all galaxies observed 
towards the cluster. Note that we will remove 'interlopers', 
non cluster members near the turn-around radius which are 
not eliminated by the 'gapper' criterion, when we calculate 
line-of-sight velocity dispersions in Section 7.1. 

For our model clusters, with their limited simulation 
volume, we did not find any gaps of 1000 km s _1 , which 
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Figure 8. Richness distribution for the ACO catalogue (thin 
line), its z < 0.08 subsample (thick line), and the rescaled 
ENACS sample (dotted line), also for z < 0.08. Numbers at the 
top indicate the traditional Abell richness classes. 

is the value used by Katgert et al. (1996) for the ENACS 
sample. So we indeed do model just what an observer would 
call the 'main' system, and we can simply apply AbelPs 
(1958) criteria to each of our simulations, without making 
any corrections. In fact, this means that we can only obtain 
the 3D richness of a cluster. 

We define the 3D richness measure C3D as the number 
of galaxies within the magnitude interval [1713, m.3 + 2] within 
a cylinder of radius 1.5/i _1 Mpc (Abell's compactness cri- 
terion) and depth 10/i -1 Mpc, which is the typical cluster 
turn-around radius for our catalogue. We chose a constant 
depth because the radius of the simulation sphere was ad- 
justed according to the expected final cluster mass in order 
to include into the modelling a sufficient volume around the 
cluster turn-around radius. 

Mazure et al. (1996) found that on average Caco is 
equal to Cro, but there exists a significant scatter in their 
relation, because for Caco a constant field correction term 
was, while for C3D the field correction term was estimated 
for each individual cluster. We used their data as plotted 
in their Fig. 1(c) to find that scatter. A linear fit gives 
"Caco-C 3 d( (:7 3d) = 0.19C 3D . 

5.2 Completeness of the model catalogue 

The main purpose of our catalogue of galaxy cluster models 
is the comparison of its statistical properties to those of the 
ENACS catalogue, so we need both samples to be complete 
in the same (observable) quantity. The ENACS sample is 
complete in the ACO richness measure CacOj so we desire 
that the model catalogue is complete in Caco as well. How- 
ever, we can only obtain Cffi for our models, so we can only 
strive for completeness in Csd- However, as Caco is statis- 
tically equivalent to C3D (Mazure et al. 1996), we can still 
compare the statistical properties of both catalogues. 

Unfortunately, completeness in C3D is not readily 
achieved. As discussed in Section 3, we had to construct 
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the catalogue using the final mass estimator Mq. Only if 
there exists a noiseless relationship between Mq and C3D 
will completeness in 3D richness immediately be guaranteed. 
However, richness is determined from the projected galaxy 
distribution, while Mq is related to the total matter distri- 
bution within a sphere. For this reason, and others which 
are discussed in Section 5.3 below, there is a rather large 
scatter in the correlation between Mq and C^o, so near 
the mass threshold some clusters drop out of the richness 
selected catalogue, while others can enter from below the 
expected mass threshold. In other words, low-mass clusters 
with a relatively high richness are, by construction, not in- 
cluded in our catalogue, and a catalogue that is complete 
in mass will not be complete in richness. As this comprises 
a significant fraction of all clusters, we need to correct for 
this. 

5.3 Correcting for incompleteness 

We correct for incompleteness in richness directly from the 
relation between richness C3D and predicted final mass 
Mq(v,x) (the catalogue defining quantity). The method 
consists of fitting an appropriate function to the {Mq^C^d) 
pairs obtained from the simulations, calculate the scatter of 
the residuals for several bins in Mq, and fit a function to 
those as well. In order to find a richness for the Mq that 
we did not select for simulation, we sample these richnesses 
from a normal distributed with mean and variance given by 
the two fits. 

In Fig. 9(a) we show the relation between C3D and Mq. 
The filled symbols indicate the simulated clusters at ag = 
0.46. The vertical dashed line shows the Mq threshold that 
selects 86 clusters, the number appropriate to the chosen 
volume (see Section 4.3), and the vertical solid line the Mq 
threshold that selects the 99 cluster models that we have 
actually simulated. The horizontal line shows the richness 
limit used to compile observational cluster catalogues (on 
average, C3D > 50 gives the same number of clusters as 
Caco > 50). 

To estimate the distribution of rich Abell clusters with 
masses below the Mq (v, x) threshold, we assume a func- 
tional form for the mass-richness relation, and measure its 
total scatter. This functional form should follow the general 
trend of the points, and stay positive for all Mq . We find 
that an exponential of the form C3D = a e bMa achieves all 
this. A fit of this form is shown in Fig. 9(a) as a solid line. 
We found a = 5.6 and b — 1.12. The scatter around this 
fit as a function of Mq is shown in Fig. 9(b), along with 
an exponential fit to it: cr(C 3D ) = 0.96 e 1ASMc . This fit 
provides the la deviations from the Mq-C^d relation that 
are plotted as dotted lines in Fig. 9(a). Combining both fits 
we derive that a{C^Y>) = 0.35C3D- 

We now assume that low-mass/high-richness clusters 
are still corresponding to initial density peaks, and should 
therefore be found amongst peaks below the Mq threshold. 
For the standard CDM scenario chosen here, the adopted 
volume of 10 7 /i~ 3 Mpc is expected to contain A r pca i ? (> 
0,^g) = 1164 peaks (see also Section 3.3.1). We then use 
both fits to estimate C3D for each Mq of the remaining 1078 
peaks below the Mq threshold. This is done by sampling 
them from a normal distribution with mean 5.6 e 112Mc 
and variance O.35C3D, he. the fits from Fig. 9. The result- 
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Figure 9. (a) 3D richness C3J3 versus predicted final cluster 
mass Mq for our simulated clusters at <rg = 0.46 (filled symbols), 
and peaks sampled below the Mq threshold for 99 clusters (open 
and grey symbols). The dashed line indicates the Mq threshold 
for 86 clusters. The horizontal lines indicate the C3D = 50 and 
C313 = 75 richness limits. The solid curve is an exponential fit 
to the points above the Mq -threshold (the vertical solid line), 
with the dotted lines showing the la deviation around that fit, 
as obtained from (b). (b) Variance in richness as a function of 
3D richness, as obtained from the simulated models, i.e. the filled 
symbols displayed in (a). The solid curve is again an expontial 
fit to these points, and was used to plot the dotted lines in (a). 

ing data points are plotted as grey dots and open circles 
in Fig. 9(a). All points below the Mq threshold that have 
C3D > 50 should be the clusters that we miss in our cata- 
logue; these points are plotted in grey. Note that most of 
them are near the Mq threshold, so it did make sense to 
run 13 more cluster models below that threshold. In fact, in 
turns out that four of those are showing up as rich clusters 
(i.e. C3D > 50) at a$ — 0.46, as can be seen in Fig. 9(a). 
According to the exponential fit, and its scatter, this is a 
plausible number. 

The distribution of all points plotted in Fig. 9(a) is the 
complete richness distribution (in a statistical sense). Since, 
on average, C3D is equal to Caco> we found the complete 
distribution function for Caco as well. Note that the num- 
ber of C3D > 50 rich Abell clusters we find from this method 
is not necessarily the amount we need (i.e. 86 in our case). 
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This will depend on the value of as, and we actually use this 
to constrain its value in Section 6.1.2. 

For the choice of a$ = 0.46 that we used to illustrate the 
method, we find that we have 78 clusters with C3D > 50, of 
which 60 were actually simulated. For the other 18 clusters 
we just sampled the richness from an estimated probability 
distribution, i.e. we did not simulate them, and therefore any 
other properties of these models are unknown. In Section 
7.1.1 we use a similar method to obtain line-of-sight velocity 
dispersions. 

5.4 A complete subsample 

If we increase the richness threshold used to define a sam- 
ple of galaxy clusters, we can obviously make the sample 
complete since the richest clusters are all simulated. Indeed 
Mazure et al. (1996) found that they could construct a truly 
complete subsample of ACO clusters for a C3D richness limit 
of 75, in a cone out to z = 0.1. After inspection of Fig. 9(a) 
we find that our simulation set at a"8=0.46 is also complete 
for C3D > 75 (i.e. above the horizontal dotted line in the 
Figure) . This allows a direct comparison of statistical prop- 
erties of these complete observational and model samples, 
but the sample size is obviously rather small: 33 and 29 
clusters respectively (for slightly different volumes, and for 
ct 8 = 0.46). 

5.5 Sources of incompleteness 

In this section we estimate the intrinsic noise of the richness 
measure due to its very definition and due to projection ef- 
fects. This allows us to find the scatter in the Mq-C^d 
relation, as the total scatter found in Section 5.3 is the com- 
bination of all these. 

5. 5. 1 Evolution of the third brightest galaxy 

Recall that richness is defined as the number of galaxies 
within the magnitude range [7713 , 7713 + 2] . Because the lumi- 
nosity function is not very well sampled at the high-end, 7713 
can vary quite significantly from cluster to cluster for equal 
total galaxy numbers and identical luminosity functions. In 
other words, the upper limit of the magnitude range, the 
magnitude of the third brightest galaxy, is far more impor- 
tant for the richness measure than the lower limit. A con- 
sequence of this is that the richness of a cluster can change 
significantly with time if 7773 changes due to, for example, a 
galaxy merger event. 

We estimate the r.m.s. scatter due to the evolution of 
7773 by taking for i = {1,2,3,4,5} and calculating the 
standard deviation of the five corresponding richness mea- 
sures. This gives us a measure for the 'jump' in richness 
due to a change of the third brightest galaxy. We find a 
spread of around 18 per cent on average, with a minimum 
of 1 per cent (model 4), and a maximum of 52 per cent 
(model 21). More quantitatively, the scatter due to 7773 as 
a function of C3D is plotted in Fig. 10(a), and a linear fit 
gives a m3 (C 3D ) = 2 + 0.18C 3D . 

5.5.2 Projection effects 

Since richness is a 'projected' quantity, measured within 
a cylinder, it is strongly dependent on the orientation of 
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Figure 10. Noise in the 3D richness measure, (a) Scatter in 3D 
richness due to 7713. (b) Scatter in 3D richness due to projection. 

the simulation volume towards the observer. This projec- 
tion effect is strongest for hierarchical formation scenarios 
like CDM, which produce highly irregular, clumpy galaxy 
distributions. We calculated the variation of richness due 
to orientation by observing each model cluster from 200 
random directions. We find that the spread in R due to 
projection effects is on average 13 per cent, with a max- 
imum of 39 per cent for model 97 (a poor cluster), and 
a minimum of 7 per cent for model 76 (a fairly compact, 
C3D=75 cluster). We find from Fig. 10(b) that, on average, 
<Voj(C* 3D ) = 5.7 + 0.072C 3D . 

5.5.3 Scatter in the Mq-C^d relation 

The noisyness in 7773 and the strong dependence on projec- 
tion make richness a measure with a rather large intrinsic 
spread. Quadratically adding these two noise sources, we es- 
timate the total observational scatter in C3D as a[ nt (C^) = 
6 + 0.2C*3D. For C"aco we need to add the scatter in 
the CACO-C3D relation as well, yielding a int (C A co) = 
6 + 0.27Caco 

Note that there are further sources of noise in richness. 
For example, one might easily miss a few galaxies when ob- 
serving a rich cluster in projection due to overcrowding, es- 
pecially if there is a large cD-galaxy present. However, we 
consider these noise sources as insignifant as compared to 
those discussed, as only 25 per cent of all clusters has a cD 



A catalogue of cluster models 15 



galaxy (Leir & van den Bergh 1977). 

With <t(C3d) known from the fit shown in Fig. 9(b), 
we can quadratically subtract <7i n t(C3D) from that to find a 
scatter of about 0.25C*3D in the Mq-C^d relation. 

6 THE FINAL CLUSTER CATALOGUE 
6.1 Timing of the models 

The normalisation of the standard Qq = 1 CDM density 
fluctuation is a free parameter that can be changed with- 
out changing the characteristic scales of the CDM scenario. 
So a numerical simulation adopting such a scenario can be 
run to an arbitrary time and then renormalised to fix the 
present epoch. Usually one assumes a linear bias between 
the galaxy and the dark matter fluctuation spectra, and the 
observation that at present the r.m.s. of the 8/i _1 Mpc Top- 
Hat smoothed galaxy distribution is about unity (Davis & 
Peebles 1983). A major advantage of actually forming galax- 
ies during the simulation is that it allows us to make a di- 
rect comparison to observations without having to assume 
a value for the bias. 

Initial conditions for the cluster models were set up such 
that as, the r.m.s. of the 8/i -1 Mpc Top-Hat smoothed total 
density field, is unity at the present epoch, i.e. there is no 
linear bias. For the standard CDM scenario as was found to 
be significantly smaller than unity in most earlier work (e.g. 
Davis et al. 1985; Frenk et al. 1990; Bertschinger & Gelb 
1991). We ran all models up to as = 0.65, and a subset to 
as — 1, the unbiased normalisation. 

6.1.1 Galaxy autocorrelation function 

A traditional measure for timing numerical simulations of 
the large-scale structure of the universe is the galaxy auto- 
correlation function. Interestingly, this function is claimed 
to be mainly determined by galaxies within 10 ft _1 Mpc of 
cluster centres (Peebles 1974; McGill 1991). The galaxy 
autocorrelation function obtained from our cluster models 
might therefore be representative of the global one. How- 
ever, quite a few problems might spoil this idea in practice. 
First of all, the dynamical evolution of both the galaxy and 
the dark matter component inside the cluster is simulated 
only approximately, because several physical processes are 
modelled only crudely, or not at all. Secondly, we did not 
transform possible cD galaxies into single galaxy particles. 
Since about 25 per cent of all observed clusters contain one 
or more cD's (Leir & van den Bergh 1977), we need to recog- 
nise these cases and separate cluster dark matter particles 
and those belonging to the cD. 

We therefore believe a determination of the galaxy au- 
tocorrelation function from our cluster model will not be 
sufficiently unbiased. For this reason we decided to use the 
normalisation found by van Kampen (1996) from several 
unconstrained simulations, i.e. average patches of universe 
for the same scenario. These 'field models' provide better 
estimates for the large-scale autocorrelation function corre- 
sponding to the cosmology chosen. The same galaxy for- 
mation parameters were used as for the cluster models pre- 
sented in this paper. The galaxy autocorrelation function 
was calculated for the simulated galaxy distribution, and 
compared to the observed one. It was found that both its 



slope and the amplitude argue for a as in the range 0.46 to 
0.56 (van Kampen 1996), although its slow evolution does 
not rule out values somewhat outside of this range. 

Once we have chosen as, i.e. defined the present epoch 
in the simulations, we can test statistical properties of the 
catalogue. Alternatively, we can use such statistics as an ad- 
ditional and independent means of normalising the cosmo- 
logical power spectrum. Besides the autocorrelation func- 
tion we use the distribution of richness as an independent 
measure to constrain as by matching it to that observed 
for Abell clusters. If we have selected the correct scenario 
and if our modelling is sufficiently realistic, both measures 
should yield consistent timing, certainly if these measures 
are on roughly the same length-scale. Note that the results 
of COBE demand as ~ 1.3 for this scenario (Bunn, Scott 
& White 1995), i.e. a linear anil-bias, but remember that 
COBE measures power on a scale well outside the dynam- 
ical range of our simulations (indicated by the two vertical 
lines in Fig. 2). This will be discussed further in Section 8. 

6.1.2 Richness distribution 

We now combine all richness determinations for the individ- 
ual models to obtain a statistical distribution function. Fig. 
11 shows histograms of the richness measure C3D found from 
our models (thin solid lines) for four values of as- Only the 
86 clusters above the Mq threshold are used for each his- 
togram. We applied the method of Section 5.3 to correct for 
incompleteness. The corrected distributions are plotted in 
Fig. 11 as thick solid lines. 

The dotted lines denote the distributions of richness 
measure Caco f° r the ENACS sample, which is complete 
(see Section 5.1). To guide the eye, we also included, as a 
dashed line, the fit to the richness distribution for the com- 
plete ACO catalogue (see Fig. 8), rescaled to the ENACS 
volume. As the distributions for Caco an d C3D are sta- 
tistically equivalent, we can directly compare our corrected 
C3D distribution to Caco distribution of the ENACS sam- 
ple and the rescaled ACO catalogue. 

We clearly observe a monotonous evolution of our cor- 
rected C3D distribution. Those for as = 0.46 or as = 0.51 
fit the observed distribution best, although the values of 0.41 
and 0.56 are not strongly ruled out. However, even smaller 
or higher values for as are excluded (we did not plot these 
in Fig. 11). We can summarize the evolution by plotting 
the density of rich clusters as a function of as- We do this 
for both the corrected C3D > 50 sample and the complete 
C3D > 75 subsample. These are shown in Fig. 12, along 
with exponential fits to them. This functional form proba- 
bly fits best for the range of as we considered here, although 
one expects the curve to stop rising exponentially for large 
as- 

From the exponential fits we find that as = 0.46 
matches best. White, Efstathiou & Frenk (1993) find a 
value of ~ 0.52 from a comparison of cluster abundances 
calculated using the Press- Schechter (1974) formalism and 
standard N-body modelling (i.e. no modelling of galaxy for- 
mation) to X-ray data. A similar comparison to the optical 
cluster abundance obtained from the median velocity dis- 
persion of Caco > 50 clusters gives erg « 0.62. Eke, Cole & 
Frenk (1996) find 0.50T0.04 from similar simulations (but 
at a higher resolution) as White et al. (ibid.), when com- 
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Figure 11. Uncorrected distributions of the richness measure 
C3J3 obtained for our model clusters (thin solid lines), and those 
corrected for incompleteness (thick solid lines), for four different 
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Figure 12. Number density of rich Abcll clusters as a function 
of erg for the standard CDM scenario. Closed symbols are for 
the corrected C3J3 > 50 sample, open ones are for the complete 
C3D > 75 subsamplc. The horizontal dashed lines indicate the 
observed densities for both samples, whereas the curves show ex- 
ponential fits to the plotted points. The preferred value for erg is 

0. 46. 

pared to just the X-ray data, where they used the same 
X-ray data as White et al. but with a slightly different way 
of inferring the X-ray cluster abundance. Unfortunately it 
is not straightforward to compare these inferred values for 
erg to the one we infer from the richness distribution, as 
there might be systematic differences between the richness 
distribution and the X-ray temperature function. 

6.1.3 Choice for erg 

The galaxy autocorrelation function and the richness distri- 
bution for the models leave a reasonably narrow range of 
allowed values for erg. In denning the catalogue using the 
expected final mass Mq, we have adopted erg = 0.46 as the 
final epoch, as it is the preferred value for the richness dis- 
tribution, which is the most restrictive of the two measures. 

By adopting erg = 0.46 we have a single catalogue that 
can be used for subsequent study. This should include a 
study of the distribution of intrinsic properties of clusters, 
and a more a detailed comparison to the ENACS and other 
observed samples. 

6.2 Selection of catalogue entries 

Now that we 'calibrated' erg to be 0.46, we select those clus- 
ter models from the set of 99 simulated that have C3D > 50 
in order to construct a richness selected catalogue. We find 
that we have simulated 60 rich clusters out of the 86 desired, 

1. e. we are 70 per cent complete in richness 630- However, 
the density of rich clusters we find for erg = 0.46 is only 
7.8 x 10~ 6 /i 3 Mpc -3 , while that should be 1.1 times larger 
according to the fit of Fig. 12 and the density found for the 
ENACS sample (Mazure et al. 1996). So we expect to find 
66 rich clusters in our Mq selected catalogue, i.e. we expect 
to be 77 per cent complete, and to miss 20 rich clusters that 
will form from peaks below the Mq threshold for the whole 
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set of 99 models. The fraction of peaks that form a rich 
cluster at the present epoch is quite small there, so it is not 
straightforward to model those clusters (see also Fig. 9(a)). 

Fortunately, as discussed in Section 5.3, we do have a 
complete subsample of C3D > 75 clusters in our simula- 
tion set, i.e. a sample for which no statistical corrections 
are needed and for which we have each entry numerically 
modelled. 

The resulting catalogue is presented in Table 1, which 
lists the defining constraints as well as some global physical 
properties (discussed below). The cluster models that are 
not rich enough are listed in italics, and those that are in the 
complete C3D > 75 are shown in bold face. Note that the 
selection of models that are part of the catalogue is specific 
for the choice of the present epoch. An extreme example 
of this is model no. 60. Its mass within 6/i -1 Mpc almost 
doubles from as = 0.46 to ag — 0.51 because another clump 
enters the Top-Hat window. However, it is also the least rich 
model of all (because its line-of-sight is perpendicular to the 
merging direction), and therefore well below the richness 
limit. 

6.3 Global properties of the cluster models 

The properties of all models are listed in Table 1. Columns 
2 to 5 give the peak parameters used to define the clus- 
ter, which also determine the expected final cluster mass 
Mq(v, x) listed in column 6. The mass within 6/i~ 1 Mpc 
found in the simulations, used to fit (21), is given in column 
7. The Abell mass Ma, i.e. the mass within a sphere of 
1.5/!,- 1 Mpc, can be obtained for most observed clusters (al- 
though usually only approximately), and is therefore listed 
in column 8. The 3D richness measure C3D, which is on 
average equal to the ACO richness measure that usually de- 
fines an observational catalogue, is listed in column 9. Since 
its value can be quite dependent on the chosen line-of-sight, 
especially for low richness clusters, we show in column 10 
the same quantity averaged over 200 lines of sight. 

The next four columns list various velocity dispersions, 
starting with the intrinsic velocity dispersions a^ m and cr ga i, 
for the dark matter and the galaxies respectively, which are 
obtained within the 3D Abell radius and divided by v3 in or- 
der to compare them to the line-of-sight velocity dispersion. 
Column 12 gives the line-of-sight velocity dispersion o\ s 
obtained from the projected galaxy distribution. Column 13 
provides, for comparison, the same quantity averaged over 
200 lines of sight. The remaining three columns provide a 
few more useful cluster properties: the masses of the third 
and the tenth ranked galaxies, and the cluster turn-around 
radius. 



7 SOME TESTABLE PROPERTIES OF 
THE n = 1 CDM CATALOGUE 

7.1 Cumulative distribution of line-of-sight 
velocity dispersions 

The velocity dispersion of a galaxy cluster generally in- 
creases with time for most cosmological scenarios. This 
makes it a useful quantity for testing the timing of our 
CDM Q — 1 model catalogue. We compare to the ENACS 
data (Mazure et al. 1996) and observations by Collins et 
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Figure 13. (a) Cumulative distributions of the line-of-sight 
velocity dispersion within 1.0 h~ ^Mpc, normalised to the 
average cluster density. The dashed line denotes the observed 
distribution for the ENACS sample. Numbers denote the value 
of <T8 for each model curve, (b) same as (a), but for dispersion 
within 1.5 /i _1 Mpc. The non-solid lincstyles represent the 
distribution for various observed samples, (c) same as (a), but 
for the complete subsample of C3J3 > 75 clusters. 

al. (1995), taken in the context of the Edinburgh-Milano 
Cluster Redshift Survey (EDSGC for short). Both these 
samples are relatively complete and should provide the best 
comparison, especially the ENACS sample because we used 
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Table 1. Parameters and global properties of the f2g = 1 CDM model cluster catalogue for crg=0.46, with Hq=50 km s 1 Mpc 
Please refer to Section 6.3 for a detailed description of each column. 
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2 
3 
4 

5 
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8 
9 
10 
11 
12 
13 
14 
15 
16 
17 
18 
19 
20 
21 
22 
23 
24 
25 
26 
27 
28 
29 
30 
3J 
32 
33 
34 
35 
36 
37 
3S 
39 
40 

41 

42 
43 

44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 



3.957 
2.986 
2.781 
3.842 
3.306 
2.874 
3.517 
3.600 
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3.470 
3.528 
3.601 
4.507 
3.206 
3.675 
3.284 
2.560 
3.653 
3.218 
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3.004 
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9.53 
5.49 
7.23 
3.87 
6.40 
4.76 
4.03 
4.74 
3.61 
3.86 
4.56 
4.20 
4.66 
5.93 
3.58 
5.41 
7.63 
4.80 
9.58 
8.69 
6.36 
6.48 
6.21 
6.68 
6.34 
3.63 
5.60 
6.97 
3.15 
6.22 
6.17 
4.54 
4.86 
8.89 
4.38 
7.66 
6.39 
5.14 
5.17 
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5.10 
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7.48 
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5.58 
5.77 
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5.83 



5.60 
3.61 
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4.73 
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4.57 
3.45 
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2.69 
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2.93 
2.24 
3.44 
4.30 
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3.45 
5.24 
3.22 
4.63 
4.05 
2.94 
3.32 
2.93 
3.63 
2.51 
4.49 
4.57 
2.49 
5.20 
2.40 
3.04 
3.07 
4.14 
3.82 
3.25 
2.59 
4.02 
2.87 



12.73 
9.32 
8.61 
11.38 
10.31 
7.89 
11.39 
11.95 
10.97 
10.33 
11.25 
10.03 
12.91 
9.94 
11.14 
10.94 
8.67 
11.59 
10.03 
9.85 
9.00 
6.99 
11.64 
9.81 
7.96 
9.79 
10.18 
11.03 
9.26 
11.16 
8.67 
11.93 
9.24 
8.85 
6.54 
8.89 
9.94 
9.44 
8.81 
8.24 
12.91 
10.38 
11.03 
8.69 
8.95 
10.27 
9.42 
11.14 
8.74 
9.00 
10.68 
8.93 
9.52 
9.43 
10.42 



notes: (a) bold face: C 3 d > 75, italics: C 3 rj < 50. 
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Table 1. continued 
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exactly the same data analysis procedures. We also consider 
the samples of Zabludoff, Huchra & Geller (1990, ZHG from 
here on) and Girardi et al. (1993, GBGMM from here on), 
although both of these are not drawn from a complete sam- 
ple. The ZHG sample might be complete but is rather small 
(25 clusters), which makes completeness hard to establish, 
while the GBGMM sample is biased towards richer clusters 
(although an attempt has been made to correct for this). 

Clearly we have to use the line-of-sight velocity disper- 
sion, denoted by u\ o s , for a comparison to observations. 
There are quite a few methods for obtaining line-of-sight ve- 
locity dispersions from observational data, all dealing with 



the removal of foreground and background galaxies, and 'in- 
terlopers', unbound galaxies near the cluster turn-around 
radius. With these galaxies removed, the velocity disper- 
sions found are generally smaller than obtained from the 
contaminated velocity distribution (e.g. Lucey 1983; Frenk 
ct al. 1990; den Hartog & Katgert 1996; Mazure et al. 1996). 

However, as discussed in Section 5.1.1, in our models 
we do not have any foreground or background galaxies due 
to the limited simulation volume. With respect to remov- 
ing interlopers, for the most conservative schemes, like '3<r- 
clipping', we find very few interlopers in our models. For 
more strict (but physically motivated) interloper removers, 
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like those proposed by den Hartog & Katgert (1996), we in- 
deed find quite a few more, and more importantly, in roughly 
the same numbers as for observed clusters. Furthermore, the 
interlopers that are removed are indeed outside the turn- 
around radius, justifying the method of den Hartog & Kat- 
gert (ibid.). 

Because EDSGC, ZHG and GBGMM use fairly conser- 
vative interloper removers, we do not need to remove any 
galaxies from our models for a comparison to these samples. 
For a match to the ENACS clusters, we apply the interloper 
remover of den Hartog & Katgert (ibid.), and also use the 
smaller aperture of 1.0/i _1 Mpc instead of the Abell radius 
used for the other samples. This smaller aperture results 
in slightly larger line-of-sight velocity dispersions, while the 
interloper remover used for the ENACS sample produces 
significantly smaller dispersions. 

7.1.1 Correction for incompleteness 

Before we can compare our line-of-sight velocity dispersions 
to those observed, we again need to correct for incomplete- 
ness in C3D ■ We use a similar method as that used to correct 
the richness distributions (see Section 5.3), but now using 
the o~\ s ,-Mq relation. We fit an exponential function to 
the (Mq,cti o s ) pairs obtained from the simulations, and 
calculate the scatter of the residuals for several bins in Mq. 
We then sample the velocity dispersions from a normal dis- 
tributed with mean and variance given by these two fits. 
As we use the same Mq's as in the correction procedure 
for richness C3D, we can link them and obtain (C3d,ci. .s.) 
pairs. Adding these to the ones taken directly from the sim- 
ulations gives us a complete distribution over line-of-sight 
velocity dispersions. 

7.1.2 Comparison to the ENACS catalogue 

The corrected 0i. o . B . distributions are plotted as solid lines 
in Fig. 13(a) for several (78 's, along with the uncorrected 
ones (dotted lines), and the observed distribution found for 
the ENACS sample (dashed line) . We find that the observed 
distribution matches the modelled one for ag w 0.35 — 0.45, 
but has a slightly different slope. 

The second comparison we make is that of the complete 
subsample of C3D > 75 clusters to a similar complete sub- 
sample of the ENACS sample. Note that we do not have to 
apply any incompleteness corrections here, as both sample 
are complete for this richness limit. The model distributions 
for ag « 0.4 matches the observed one best, as shown in Fig. 
13(b). The shapes of the curves compare better than for the 
C3D > 50 sample. 

7.1.3 Comparison to other observed samples 

We can do the same comparison for the EDSGC, GBGMM, 
and ZHG samples, although only the first of these is really 
complete in richness. The line-of-sight velocity dispersions 
within the Abell radius for several values for 0% are plotted 
in Fig. 13(c), along with the observed distributions. Please 
note that the interloper removal method of den Hartog & 
Katgert (1996) was not used here, as it was not used for any 
of these observed distributions either, and also note that the 
Abell aperture was adopted. 



For the EDSGC sample we find that the model curve 
fits best for ag « 0.4. The other two samples prefer much 
higher values, but are not complete and biased to the richest 
clusters, and therefore to higher velocity dispersions. 

7.2 Other observables 

To test the chosen scenario one can use many other ob- 
servables besides the one considered above. For example, 
the line-of-sight velocity dispersion profiles might prove use- 
ful, although unfortunately at present these can only be ob- 
tained for a few bins. However, a classification derived from 
this binned profile can be very instructive (den Hartog & 
Katgert 1996). 

Another distribution function is that of the projected 
shapes of clusters. In a preliminary comparison presented 
elsewhere (de Theije, Katgert & van Kampen 1995) it was 
found that the predicted distribution of shape parameters 
of model clusters differs significantly from the observed one. 
Although the observed sample used in the study of de Theije 
et al. (ibid.) is not complete, the deviations are unlikely 
to be due to incompleteness. For a subset of our model 
catalogue that was rerun for the = 0.2 CDM scenario (in 
future work we will build complete catalogues for this and 
other scenarios) there seems to be better agreement between 
models and data. 



8 SUMMARY AND DISCUSSION 

The aim of this paper was to present a technique that 
produces a fair catalogue of individually simulated high- 
resolution cluster models. We used this technique to build 
one specific model catalogue that we compared to observed 
catalogues. 

Under the assumption that rich clusters of galaxies 
originate from peaks in the initial density field Gaussian 
smoothed at 4/i _1 Mpc, we select a catalogue on the basis 
of characteristic parameters of these peaks. Given the fact 
that we need to use the constrained random field method 
to generate initial conditions for the individual numerical 
models, we need to select the catalogue on the basis of lin- 
ear functionals of the initial density field. With this in mind, 
we argued that cluster mass is the best catalogue denning 
quantity, because it can be predicted reasonably well from 
a linear function of the amplitude and the curvature of the 
initial density peak which is the progenitor of the cluster. 
Besides this practical argument, total cluster mass is a ba- 
sic property of a galaxy cluster. 

For the individual models a numerical code was used 
that contains a prescription to form galaxies, and also al- 
lows galaxies already formed to grow (i.e. accrete particles) 
and merge with other galaxies. Not only does this allow us 
to directly compare to optical catalogues of galaxy clusters, 
but the properties of the galaxy population also directly set 
the present epoch for the cosmological scenario in which the 
model catalogue is embedded, i.e. they set the amplitude of 
the initial cosmological fluctuation spectrum. Not all cos- 
mological scenarios allow this, but the scenario we chose, 
viz. fio = 1 CDM, does permit this a posteriori normali- 
sation and also makes it possible to compare our results to 
published numerical work. 
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We used a traditional, single large-scale simulation in 
order to test several issues related to the construction of the 
catalogue and the reliability of the individual models. We 
traced peaks in these simulations, and found that they do 
not merge during the entire run. So we indeed have a one-to- 
one mapping of initial peaks to final rich Abell clusters. Note 
that this may not be so in other cosmological scenarios. We 
also found that both the distribution and the evolution of 
the peak parameters in the single large-scale low-resolution 
simulation are similar to those of the set of individual high- 
resolution simulations. Finally we found that the expected 
mass relation, a linear function of peak amplitude and peak 
curvature, is almost equal to that obtained from the clusters 
extracted from the single simulation. 

Having built a catalogue selected on expected final clus- 
ter mass, we found that is was 70 per cent complete for 
C*3D > 50. This incompleteness mostly concerns the low- 
mass end and is mainly due to the noise still present in the 
relation between final cluster mass and initial peak parame- 
ters. But note that the richness measure itself has a rather 
large intrinsic noise originating from its very definition. To 
allow a comparison of statistical distribution functions of 
model clusters properties to those observed we devised a 
method to corrected for this incompleteness. However, we 
do have a complete subsample for C3D > 75, which permits 
a fully unbiased comparison to observations. 

We use the richness distribution of both samples to set 
the present epoch in the models, i.e. fix erg. A value of 0.46 
for as matches best; this is also in the range of allowed values 
derived from a match of the galaxy autocorrelation function 
for a set of field models similar to the clusters models dis- 
cussed here to the observed autocorrelation function (van 
Kampen 1996). The elimination of the free parameter o-g 
of the = 1 CDM scenario allows us to test it with other 
measures. 

We find that the cumulative distribution of line-of-sight 
velocity dispersions obtained for the ENACS sample (Kat- 
gert et al. 1996), and for other samples as well, fits the model 
distribution best for ct 8 w 0.4. This is fairly consistent with 
both the auto-correlation function and the richness distri- 
bution. We can obtain a slightly larger og if we invoke a 
velocity bias, but this is not seen in our models if we com- 
pare the intrinsic velocity dispersions of the dark matter to 
those of the galaxies (van Kampen 1994). Other studies, 
like those of Cen & Ostriker (1992) and Katz et al. (1992) 
do not find a significant velocity bias either. 

However, there are also problems with the standard 
CDM scenario. We mention two major ones here (see Os- 
triker 1993 for a more extensive discussion). De Theije et al. 
(1995) compared the distribution of cluster shapes for our 
catalogue to an observed sample, and found that these differ 
significantly. The results from COBE imply that og « 1.3 
for our chosen scenario (Bunn et al. 1995). Note however 
that COBE traces length-scales well outside the dynamical 
range of the individual models. We will need to investigate 
whether these inconsistencies remain for other tests before 
we completely discard the flo = 1 CDM scenario on the 
scales we studied here. It is still useful to test more of the 
statistical properties of this scenario, and use successes and 
failures of such tests as a guide for a new choice for the 
most likely cosmological scenario. Of the scenarios plotted 
in Fig. 2, we should take a closer look at the k~ 2 spectrum, 



or low-f2o CDM, as both of these have more power on larger 
scales and therefore are in better agreement with the COBE 
results. 
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APPENDIX A: EXPRESSIONS FROM THE 
THEORY OF GAUSSIAN RANDOM FIELDS 

In this Appendix we summarize some functions taken from 
Bardeen et al. (1986; BBKS from here on), slightly rewritten 
for our purposes, as used in Section 2. The first function is 
a factor in the probability distribution for the peak height 
v, which depends on the spectral parameters 7 and Rq. It 
depends on 7 and v only: 



7 V 



^ l(7 '^ l + C 2 (7)e-C 3 (7V 

(from equation (4.4) of BBKS), with 
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The second function is used for the curvature distribution 
function: 
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The next two functions are used for the asphericity proba- 
bilities: 
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Finally we list the function used to describe the asphericity 
of peaks in the second order Taylor expansion, which is a 
function of the spherical coordinate angles 9 and (j> and axial 
ratios ayi and 013: 

^"5(012, ai3, 0, 
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1) + (2af 2 - 1) cos z 4>} 

This is an adaptation from eq. (7.4) of BBKS. It is equal to 
unity for spherical peaks. 

We draw from the probability distributions (10) and 
(11) by using the rejection method. The comparison func- 
tion needed is found by substituting v for ^1(7, v) in (10). 
Its integral is 

V 



(2tt) 2 RI 



2- {2 + v 2 )e- v 12 



(A8) 



The comparison function for drawing x, f c (x,fi, a), is ob- 
tained by replacing Ti(x) with x 3 , where we have defined 
fj, = -yu and a = \fl — -y 2 . The integral F c (x,fi,a) of the 
comparison function is 



- 



F c (x, fj,,a) = — a 2 (^i 2 + 2<r 2 + fii + i 2 )e 

+ ^yf ( M 2 + 3. 2 ) Erf(^) 



(A9) 
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An approximate formula for the expectation value of x for 
a given v is given by eq. (6.14) of BBKS: 

M= 3(1 - 7 2 ) + (1-216 - 0.9 7 4 )e^ 3,y2 / 8 

{X) 7 " [3(1 -7 2 )+ 0.45 + (7^/2)2] 1/2 + 7I ,/2 ( j 

The expectation values for 012 and £113 are functions of y = 
(6 + 5X 2 )- 1 / 2 : 



(012> 



(ai3> 



1 - 60y 4 



1 - 3y + 30y 4 

1 + 3y - 30y~ 
l-3y + 30y 4 



(All) 



APPENDIX B: DERIVATION OF THE 
PREDICTED MASS RELATION 

In this appendix we show how the curvature of the initial 
peak, x, comes into play in the non-linear evolution of total 
cluster mass. We start from the initial density profile, given 
by (19): 



5 G (r) 



■ ix „ . . v 



■V 2 Cg(0 • (Bl) 



o"o(l~7 2 ) 3o- (l-7 2 ) 
We neglect non- linear evolution of ao, so that 7 remains 
constant (o\ and 02 remain well below unity). We assume 
that the autocorrelation function is a power law with index 
n, which evolves as n(t) = n(tj)e(t), where e(t) is unity at 
t = tj and increases slowly. This gives us approximately 

V 2 fc(0,t)^ e(f)( ;j;;jf 1 +1) a 2 (t)V 2 fc(0,t 1 ) 

= _3^)("(W) + D 7 2 g (t)i 



(B2) 
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where we have used eq. (20) for the second part, and thus 



S G (Q,t) 
with 



v(t) - -yx(t) + rj ^ 1 x ( t ) - 



(I-7 2 ) 
e(t)[n(ti)e(t) + l] 
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<r (t) (S3) 



(54) 



ra(ti) + 1 

which monotonically increases (slowly) from its initial value 
r?(ti) = 1. The second term in eq. (-B3) evolves somewhat 
faster than the first because of the increasing slope n(t) of 
the autocorrelation function. This might seem unimportant, 
but the point is that the x terms, which grow relatively fast 
(as clearly seen in Fig. 6) , do not cancel anymore because of 
the different growth rates of the two terms in eq. (B3). 
Rearranging the terms, we find 

S G {Q,t) = _ ^ 2 ao(t) . (B5) 

We can rewrite this in terms of growth factors for v and x: 
*g(0,*) = {g v (t>i+9x(t)xi}(7 (t )a(t) . (B6) 
with V[ = v{t\), x,- L = x{t\), and 



and 

gx{t) = ^~ 

X\ 
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(58) 



The CDM spectrum smoothed on 4/i _1 Mpc has an autocor- 
relation function with an initial slope of n w —1.5 for the 
range of wavenumbers in our simulations (Efstathiou 1990). 
As an example of what the final relation will look like we 
assume that 77^0) ~ 1-2. Because 7(4/i _1 Mpc) = 0.729 for 
CDM, we get 



v\ x\ 



(B9) 



The peak amplitude u(t) just grows, but the evolution of 
the curvature x(t) is more complicated. Van Haarlem & 
van de Weygaert (1993) performed cluster simulations with 
constraints set on the same smoothing scale, and found 
that a small initial x\ results in a larger final x(to) than 
a larger initial x\. For three simulations with v x — 3 
and Xi's corresponding to the 5, 50 an 95 percentiles of 
P(x\3ao, 4/i~ 1 Mpc), a linear relation provides a good fit to 
their results: 



a; (t ) « 27-4.62:; 



(B10) 



We should note that both constants are present epoch values 
of unknown functions of time, but here we just look for 
a relation between initial and final values. This relation 
is obtained for v\ = 3 only, but since that is roughly the 
average value for our model cluster sample, a linear relation 
for a:(fo) seems to be a good approximation for all v\. 

We now assume a linear relation for v(to) as well, i.e. 
^(to) w avi + b, where again the fitting parameters are func- 
tions of to- This linear relation is a good approximation 
because we only consider a relatively small range of V[, so 
that the non-linear relation between i^(to) and v\ can locally 
be well approximated by a linear function. 

Combining both linear relations we can write the ex- 
pected final total mass as 

M G (u i ,x i ,t ) = {27TR' G f /2 p[c Q +c 1 Ui + C2X i ] . (Bll) 

The constant c\ will be larger than one, as v(t) is a mono- 
tonically growing function. On the other hand, we expect 
C2 to be negative on the basis of eq. (B10). 



APPENDIX C: RELIABILITY OF THE 
CLUSTER MODELS 

Our models could be susceptible to various systematic effects 
and errors. The first source of error that comes to mind is 
the N-body code itself, since the integration of the particle 
trajectories is performed with a certain accuracy. However, 
the energy is conserved sufficiently well for all models: about 
one procent change in E/U(0) (the total energy divided by 
the initial total potential energy). The shape of the initially 
spherical distribution of particles changes little and its bor- 
ders do not move inwards. This means that we have indeed 
included the complete region of influence of the cluster, i.e. 
the turnaround radius is well within the simulation volume. 
The autocorrelation function lacks some of the large-scale 
power at late times, but that is only the case for ag larger 
than the value of 0.46 that we adopt. However, we will need 
to take larger simulation volumes for scenarios with more 
power on large scales, like the CDM fio — 0.2 scenario. 

Systematic errors are most likely due to the approxi- 
mate way in which we model the formation and evolution 
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of galaxies. For example, mass loss of galaxies is not mod- 
elled, and dynamical friction may therefore be unrealistically 
large. However, the galaxy formation parameters were cho- 
sen such that only the cores of particle groups that were 
found to be virialised were tranformed into a single galaxy 
particle (see also van Kampen 1995). This allows the re- 
maining outer particles to be stripped, even though this need 
not necessarily happen. The crude way in which merging is 
modeled is probably the biggest defect. Merging may there- 
fore be underestimated in some regions and overestimated in 
others. In future work we will try to add a separate merger 
criterion, even though this will introduce more parameters 
to the modelling. 



